/*
 * 
 *
 * Authors: Soeren Gebbert, soeren.gebbert@vti.bund.de
 *          Axel Don, axel.don@vti.bund.de
 *
 * Copyright:
 *
 * Johann Heinrich von Thünen-Institut
 * Institut für Agrarrelevante Klimaforschung
 *
 * Phone: +49 (0)531 596 2601
 *
 * Fax:+49 (0)531 596 2699
 *
 * Mail: ak@vti.bund.de
 *
 * Bundesallee 50
 * 38116 Braunschweig
 * Germany
 *
 *
 * This program is free software; you can redistribute it and/or modify
 * it under the terms of the GNU General Public License as published by
 * the Free Software Foundation; version 2 of the License.
 *
 * This program is distributed in the hope that it will be useful,
 * but WITHOUT ANY WARRANTY; without even the implied warranty of
 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PUrunOSE.  See the
 * GNU General Public License for more details.
 */

#include <vector>
#include <iostream>
#include <stdlib.h>
#include <math.h>
#include <stdio.h>
#include <fstream>
#include <time.h>
/* The ANN - approximate nearest neighbor searching library
 * of version 01/27/10 (Version 1.1.2) is needed.
 * ANN.h and libANN.a are required
 * 
 * http://www.cs.umd.edu/~mount/ANN/
 * 
 */
#include "ANN.h"

#include "carb_particle.h"
#include "micr_cluster.h"
#include "initial_parameters.h"
#include "runtime_parameters.h"
#include "event_base.h"
#include "output.h"

using namespace std;

static int ANN_radius_search(vector<micr_cluster*> &mc, vector<carb_particle*> &cp);

/* -------------------------------------------------------------------------- */
int main(int argc, char **argv) {
    unsigned int timestep;
    unsigned int i, j;

    ofstream run_file;
    vector<carb_particle*> cp;
    vector<micr_cluster*> mc;
    vector<runtime_parameters*> runparams;

    /* Create initial parameter */
    initial_parameters* init = new initial_parameters();
    /* Read values from file if present */
    if (argc > 1)
        init->read(argv[1]);

    /* Create runtime paramater */
    runtime_parameters* runparam = new runtime_parameters();
    /* Read values from file if present */
    if (argc > 2)
        runparam->read(argv[2]);

    /* We have a runtime parameter for each timestep */
    for (i = 0; i <= init->number_of_timesteps; i++)
        runparams.push_back(runparam);

    /* Initiate the random number generator */
    if(init->seed == 0) {
        srand48(time(NULL));
    } else {
        srand48(init->seed);
    }
    
    /* initialize the micro-organism cluster */
    for (i = 0; i < init->num_micr_cluster; i++) {
        micr_cluster *c = new micr_cluster(init);
        mc.push_back(c);
    }

    /* Set up the output run file*/
    run_file.open(init->run_filename.c_str());

    /* Write header*/
    run_file << "timestep,micr_mass,in,eaten,out,carb_total_mass,carb_num" << endl;

    /* Initial valuest */
    double mass = 0.0;
    for (i = 0; i < mc.size(); i++)
        mass += mc[i]->get_mass();

    run_file << 0 << "," << mass << ",0,0,0," << cp.size() << endl;
    /* Main loop */
    for (timestep = 1; timestep <= init->number_of_timesteps; timestep++) {
        int in = 0;
        int eaten = 0;
        
        /* The active runtime parameter */
        runtime_parameters* run = runparams[timestep];

        /* Generate new carbon particles */
        for (i = 0; i < run->new_carb_particles; i++) {
            carb_particle *p = new carb_particle(init, run);
            cp.push_back(p);
            in++;
        }

        /* ********** *
         * ** MOVE ** *
         * ********** */

        /* Set the runtime parameter and move the carbon particles */
        for (i = 0; i < cp.size(); i++) {
            cp[i]->set_runtime_parameters(run);
            cp[i]->move();
        }

        /* Set the runtime parameter and move the micro organism */
        for (i = 0; i < mc.size(); i++) {
            mc[i]->set_runtime_parameters(run);
            mc[i]->move();
        }

        /* ********* *
         * * EVENT * *
         * ********* */

        /* Iterate over each event */
        for(vector<event_base*>::iterator it = init->events.begin(); it !=  init->events.end(); ++it)
            (*it)->modify_particles(timestep, init, run, mc, cp);

        /* ********* *
         * ** EAT ** *
         * ********* */
        
        if(init->search_method.find("ANN") != string::npos) {
            eaten = ANN_radius_search(mc, cp);
        } else {
            /* Search neighbour carbon particles using brute force */
            for (i = 0; i < mc.size(); i++) {
                /* We compute the radius at this point or performance reasons */
                float radius = mc[i]->get_radius();

                for (j = 0; j < cp.size(); j++) {
                    /* Ignore inactive particles */
                    if (cp[j]->is_active() == false)
                        continue;
                    /* Check radius */
                    if (mc[i]->is_particle_in_radius(cp[j], radius) == false)
                        continue;
                    /* Micro-organism eat the active carbon in range */
                    mc[i]->eat(cp[j]->get_eaten());
                    eaten++;
                }
            }
        }

        /* *************** *
         * ** RESPIRATE ** *
         * *************** */

        /* Organism respiration */
        for (i = 0; i < mc.size(); i++)
            mc[i]->respirate();


        /* ************** *
         * ** CLEAN UP ** *
         * ************** */

        /* Count age and remove un-needed carbon paricles */
        int out = 0;
        vector<carb_particle*> cp_tmp;

        for (i = 0; i < cp.size(); i++) {
            if (cp[i]->is_active()) {
                cp[i]->aging();
                cp_tmp.push_back(cp[i]);
            } else {
                out++;
                delete cp[i];
            }
        }

        double carb_mass = 0.0;
        cp.clear();
        /* Copy back the active carbon particles and compute the total mass */
        for (i = 0; i < cp_tmp.size(); i++) {
            cp.push_back(cp_tmp[i]);
            carb_mass += cp_tmp[i]->get_mass();
        }
        cp_tmp.clear();

        if (cp.size() > init->max_num_carb_particle)
            cerr << "WARNING: the maximum number of carbon particles are reached: " << cp.size() << endl;

        /* ************ *
         * ** OUTPUT ** *
         * ************ */

        double mass = 0.0;
        for (i = 0; i < mc.size(); i++)
            mass += mc[i]->get_mass();

        run_file << timestep << "," << mass << "," << in << "," << eaten
                << "," << out << "," << carb_mass << "," << cp.size() << endl;

        /* Write output files */
        if (timestep % init->output_step == 0) {
            cout << timestep << " mass " << mass << " in " << in << " eaten "
                    << eaten << " out " << out << " carb_total_mass " << carb_mass << " carb_num " << cp.size() << endl;
            write_carb_particle_output_files(timestep, init, cp);
            write_micr_cluster_output_files(timestep, init, mc);
            write_profile_output_file(timestep, init, mc, cp);
        }
    }

    run_file.close();

    delete init;
    delete runparam;

    annClose(); // done with ANN
    
    return 0;
}

/* -------------------------------------------------------------------------- */

int ANN_radius_search(vector<micr_cluster*> &mc, vector<carb_particle*> &cp) {
    ANNpointArray dataPts; // data points
    ANNpoint queryPt;      // query point
    ANNidxArray nnIdx;     // near neighbor indices
    ANNdistArray dists;    // near neighbor distances
    ANNdist sqRad;         // Squared radius
    float radius;
    ANNkd_tree* kdTree;    // search structure
    unsigned int nPts = cp.size();  // actual number of data points
    unsigned int dim = 3;
    unsigned int i, j;
    unsigned int eaten = 0;
    float x, y, z;


    queryPt = annAllocPt(dim);        // allocate query point
    dataPts = annAllocPts(nPts, dim); // allocate data points
    nnIdx = new ANNidx[nPts];         // allocate near neigh indices
    dists = new ANNdist[nPts];        // allocate near neighbor dists
    
    /* Fill the data points */
    for(i = 0; i < cp.size(); i++) {
        x = cp[i]->get_x();
        y = cp[i]->get_y();
        z = cp[i]->get_z();
        dataPts[i][0] = x;
        dataPts[i][1] = y;
        dataPts[i][2] = z;
    }

    /* Create the kd-tree */
    kdTree = new ANNkd_tree(dataPts, nPts, dim);

    /* Search carbon particles in radius of each micro-organism */
    for(i = 0; i < mc.size(); i++) {

        /* Build query point */
        x = mc[i]->get_x();
        y = mc[i]->get_y();
        z = mc[i]->get_z();
        queryPt[0] = x;
        queryPt[1] = y;
        queryPt[2] = z;
        
        /* Compute radius */
        radius = mc[i]->get_radius();
        sqRad = radius*radius;
        /* Search nearest neighbour */
        kdTree->annkFRSearch(queryPt, sqRad, nPts, nnIdx, dists, 0.0); 
        
        /* Iterate over each carbon particle found */
        for(j = 0; j < nPts; j++) {
            /* Break if the index is invalid */
            if(nnIdx[j] == ANN_NULL_IDX)
                break;
            /* Break if this particle was already eaten */
            if(cp[nnIdx[j]]->is_active() == false)
                continue;
            /* Eat ! */
            mc[i]->eat(cp[nnIdx[j]]->get_eaten());
            eaten++;
        }
    }
    
    /* Clean up */
    delete [] nnIdx;
    delete [] dists;
    annDeallocPt(queryPt);
    annDeallocPts(dataPts);
    delete kdTree;

    return eaten;
}
